Density-functional study of the structure and stability of ZnO surfaces 



O 

o 

(N 
G 

1—5 

(N 



B. Meyer and Dominik Marx 
Lehrstuhl fur Theoretische Chemie, 
Ruhr-Universitdt Bochum, 44780 Bochum, Germany 
(February 1, 2008) 

An extensive theoretical investigation of the nonpolar (1010) and (1120) surfaces as well as the polar 
zinc terminated (0001)-Zn and oxygen terminated (000l)-O surfaces of ZnO is presented. Particular 
attention is given to the convergence properties of various parameters such as basis set, k-point 
mesh, slab thickness, or relaxation constraints within LDA and PBE pseudopotential calculations 
using both plane wave and mixed basis sets. The pros and cons of different approaches to deal 
with the stability problem of the polar surfaces are discussed. Reliable results for the structural 
relaxations and the energetics of these surfaces are presented and compared to previous theoretical 
and experimental data, which are also concisely reviewed and commented. 
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I. INTRODUCTION 

The II- VI semiconductor ZnO has become a frequently 
studied material in surface science because of its wide 
range of technological applications. ZnO is a basic ma- 
terial for varistors, thyristors and optical coatings. In 
addition, its direct band gap makes it an interesting can=. 
didate for blue and UV emitting LEDs and laser diodes.El 
The electronic and structural properties of the ZnO sur- 
faces are in particular important in its applications as 
chemical sensor in gas detecting systems and as cata- 
lyst for hydrogenation and dehydrogenation reactions. In 
combination with Cu particles at the surface, ZmO is a 
very efficient catalyst for the methanol synthesiscl where 
it is employed in industrial scale. The mechanism behind 
the enhanced catalytic activity when combined with Cu 
is poorly understood. However, before this interesting 
interplay between the ZnO substrate and the Cu parti- 
cles can be addressed, a thorough understanding of the 
underlying clean ZnO surfaces is necessary. 

From a physical/chemical point of view, ZnO is a very 
interesting material because of the mixed covalent/ionic 
aspects in the chemical bonding. ZnO crystallizes in 
the hexagonal wurtzite structure (B4) which consists of 
hexagonal Zn and O planes stacked alternately along the 
c-axis (see Fig. [|). Anions and cations are 4-fold coordi- 
nated, respectively, like in the closely related zincblcndc 
structure. A tetrahedral coordinated bulk structure is 
typical for rather covalent semiconductors. On the other 
hand, ZnO shows great similarities with ionic insulators 
such as MgO.B This is why ZnO is often called the 'ionic 
extreme' of tetrahedral coordinated semiconductors. 

Wurtzite crystals are dominated by four low Miller 
index surfaces: the nonpolar (1010) and (1120) sur- 
faces and the polar zinc terminated (0001)-Zn and the 
oxygen terminated (0001)-O surfaces (see Fig. |l|). By 
ion sputtering and annealing at not too high tempera- 
tures all four surfaces can be prepared in a bulk termi- 



nated, unreconstructed state, where the surface atoms 
only undergo symmetry conserving relaxations. A typ- 
ical p(lxl) pattern is observed in low-energy electcpn. 
diffraction (LEED) and other diffraction experiments. lTEI 
Although in a recent He-scattering experiment^ it was 
shown that O-terminated (0001) surfaces with p(lxl) 
LEED patterns are usually hydrogen covered whereas the 
clean O-terminated surface exhibits a (3x1) reconstruc- 
tion, we will focus in this study on the clean, unrecon- 
structed surfaces of ZnO. 

In the present paper, we investigate all four main 
crystal terminations of ZnO. The fully relaxed geomet- 
ric structures and the surface/cleavage energies have 
been calculated using a first-principles density-functional 
(DFT) method. We have employed both, a local-density 
(LDA) and a generalized-gradient approximation (GGA) 
functional. We will discuss the relative stability of the 
four surfaces and how the surface relaxations of the non- 
polar faces are connected to the covalency /ionicity of the 
chemical bond in ZnO. Finally, a detailed comparison 
with existing theoretical and experimental results will be 
given. 

The nonpolar (1010) surface of ZnO has been the focus 
of several experimental and theoretical studies. However, 
the form of the relaxation oLthe surface atoms is still 
very controversial. Duke et al.tj concluded from their best 
LEED analysis! that the top-layer zinc ion is displaced 
downwards by Adj_(Zn)=— 0.45±0.1 A and likewise the 
top-layer oxygen by Ad±(0)=— 0.05±0.1 A, leading to a 
tilt of the Zn-0 dimer of 12°±5°. No compelling evidence 
for lateral distortions within the first layer or for second- 
layer relaxations were obtained, but small improvements 
could be achieved by assuming a lateral displacement of 
the Zn ion toward oxygen by Ad\\ (Zn)=0.1±0.2 A.El The 
strong inward relaxation of the Zn ion was later con- 
firmed by Gopel at alO in an angle-resolved photoemis- 
sion experiment. By comparing the relative position of 
a particular surface state with its theoretically predicted 
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FIG. 1. Wurtzite structure (B4) of ZnO with the po- 
lar zinc terminated (OOOl)-Zn, the polar oxygen terminated 
(0001)-O, and the nonpolar (1010) surfaces. 

geometry dependence, a Zn displacement downwards by 
Ad_L(Zn)=— 0.4 A was concluded. 

In contrast, Jedrecy et al.EHl found best agreement with 
their gracing incidence X-ray diffraction data (GIXD) for 
a structural model where the top-layer zinc atom is dis- 
placed downwards by only Aci^(Zn)=— 0.06±0.02 A and 
shifted toward oxygen by Ady (Zn)=0.05±0.02 A. How- 
ever, for their samples they observed a high density of 
steps and from their best-fit model they predict rather 
high vacancy concentrations in the first two surface lay- 
ers with occupancy factors of 0.77±0.02 and 0.90±0.04 
for the first and second layer, respectively. On the other 
hand, Parker et al£3 reported scanning tunneling mi- 
croscopy (STM) images of the nonpolar (1010) surfaces 
with atomic resolution where large flat terraces are found 
and no defects are visible in areas as large as 11x14 sur- 
face unit cells. Due to the small scattering contribution, 
the position of oxygen could not be determined very ac- 
curately in the GIXD experiment of Ref. |lj. The result 
of the best fit was that O relaxes further toward the bulk 
than Zn with Ad_L(O)=-0. 12+0.06 A. This would be 
very unusual since to our knowledge no (1010) wurtzite 
or (110) zincblende surface structure has been reported 
where the surface dimers tilts with the cation above the 
anion. 

First theoretical investigations of the (1010) surface 
were done using empirical tight-binding (TB) models^ 
With two very different TB models Wang and Duketa 
found a strong displacement _of Ad^(Zn)=— 0.57 A, 
whereas Ivanov and PollmanrJlj obtained an almost 
bulk-like surface geometry. A recent calculation with 
atomistic potentials based on a shell modelll3 predicted 
Adj_(Zn)=— 0.25 A and a rather strong upward relax- 
ation of the second-layer Zn of +0.165 A—. 

Several ab-initio studies (DFT-LDA,EZI Hartree-Fock 
(HF)Jl3 and a hybrid HF and DFT method using the 



B3LYP functional^) employing Gaussian orbitals as ba- 
sis functions to solve the electronic structure problem 
favor small inward relaxations of Zn and small tilts of 
the ZnO-dimers of 2°-5°. However, it is questionable, 
if these studies represent fully converged results. There 
is only one recentr-first-principles DFT-LDA calculation 
using plane wavesE3 where larger relaxations with a tilt 
of 11.7° were obtained. 

The nonpolar (1120) ZnO surface has been less fre- 
quently studied than its,_Ll010) counterpart. The 
two tight-binding modelai-ftLa predicted the same relax- 
ation behavior for the_(1120) as for the (1010) sur- 
face: Wang and Dukaifl found a strong zinc displace- 
ment of Adj_(Zn)=— 0.54 A toward the bulk whereas the 
TB model of Ivanov and Pollmann preferred an almost 
bulk-like surface structure. With a first-jarinciples hy- 
brid B3LYP method Wander and Harrisono found much 
smaller relaxations for the (1120) surface than for the 
(1010) face, but not all degrees of freedom were relaxed 
in this study. To our knowledge there has been no quan- 
titative experimental investigation. 

Coming to the polar surfaces, we encounter the fun- 
damental problem that in an ionic model these surfaces 
are unstable and should-4iot exist. They are so-called 
'Tusker type 3' surfaces,c3 and with simple electrostatic 
arguments it can be shown that, the surface energy di- 
verges for such a configuration.!^ To stabilize the polar 
surfaces a rearrangement of charges between the O- and 
the Zn-terminated surfaces needs to take place in which 
the Zn-terminated side becomes less positively charged 
and the O-terminated face less negative. In fact, most 
polar surfaces show massive surface reconstructions og 
exhibit facetting to accommodate the charge transfer.Q 
Also randomly distributed vacancies, impurity atoms in 
the surface layers, or the presence of charged adsorbates 
are possible mechanisms to stabilize polar surfaces. How- 
ever, the polar ZnO surfaces are remarkably stable, and 





LDA 


PBE 


Expt. 


a [A] 


3.193 (-1.7%) 


3.282 (+1.0%) 


3.250 


c[A] 


5.163 (-0.8%) 


5.291 (+1.6%) 


5.207 


c/a 


1.617 


1.612 


1.602 


u 


0.3783 


0.3792 


0.3825 


Bo [GPa] 


161 


128 


143 


Pt [GPa] 


9.0 


11.8 


9.0-9.5 



TABLE I. Computed and experimental values of the struc- 
tural parameters for bulk ZnO. a and c are the lattice con- 
stants, u is an internal coordinate of the wurtzite structure 
which determines the relative position of the anion and cation 
sublattice along the c axis, Bo is the bulk modulus, and pt is 
the transition pressure between the wurtzite (B4) and rock- 
salt (Bl) structure of ZnO. Experimental values are from 
Refs. pj[-|24l Relative deviations from experiment are given 
in parenthesis. 
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(b) Side view: 




FIG. 2. Schematic diagram of the surface geometry and 
the independent structural parameters of the nonpolar (1010) 
surface. The brackets indicate the two atomic layers shown 
in top and side view. Open and filled symbols are the O 
and Zn ions, respectively, and the solid lines represent near- 
est-neighbor bonds. The atoms in the first layer are shown 
by solid/black, second layer atoms by dashed/shaded circles. 
The surface unit cell is indicated by dashed lines. 

many experiments suggest that they are. in an unrecon- 
structed, clean and fully ordered stateE Despite many 
investigations it is still an open question how the polar 
ZnO surfaces are stabilized.El 

Assuming clean and unreconstructed surfaces, the re- 
duction in surface charge density can only occur from a 
redistribution of the electrons. Negative charge has to 
be transferred from the O-terminated face to the Zn- 
terminated side, leading to partially occupied bands at 
the surface. This so called 'metallization of the surfaGe^. 
has been used by all previous ab-initio calculationscjOEl 
to model the polar ZnO surfaces and will also be em- 
ployed in the present study. However, whether or not 
the surfaces are metallic will depend on the width of the 
partially occupied bands. From another point of view, if 
the polar surfaces were stabilized by vacancies, defects or 



adsorbates, many defect states would be created. Now 
if we think of somehow averaging over the surface, the 
defect states would form a partially filled band. In this 
sense, the 'metallization' may be regarded as some 'mean- 
field' description for a situation where many defect states 
are present. 

Several attempts have been made to determine the 
layer relaxations of the unreconstructed polar sur- 
faces. In an_ early dynamical LEED analysis Duke 
and LubinskyEj found an outer Zn-0 double-layer spac- 
ing of <ii2=0.607A for the Zn-terminated surface and 
c?i2=0.807A for the O-terminated face. Unfortunately, 
this analysis was based on an early bulk structure of ZnO, 
see Ref. |29|, in which the bulk double-layer spacing was 
assumed to be 0.807 A instead of 0.612 A. 

For the Zn-terminated surface, it was concluded from 
the comparison of X-ray ohoto-diffraction (XPD) data 
with scattering simulationsE3 that any inward relaxation 
of the surface Zn layer can be ruled out. Coaxial impact- 
collision ion-scattering spectroscopy!^! (CAICISS) pro- 
posed an expansion of d\2 by +0.35 A. Also an expan- 
sion of di2 by +0.05 A for the Zn-terminated surface was 
found in an GIXD measurement H3 In this experiment, 
the X-ray data could best be fitted by assuming a ran- 
dom removal of 1/4 of the Zn atoms in the surface layer. 
On the other hand, from the shadowing and blocking 
edges of a low-energy alkali ion scatterings (LEIS) exper- 
iment no evidence for substantial quantities of point de- 
fects in the Zn-terminated as well as in the O-terminated 
surface was found. 

Foc-the O-terminated surface, it was concluded from 
LEISo that the Zn-0 double-la-wr spacing dyi is close 
its bulk value. An XPD studya found a contraction 
of 25% of di2, but like in the LEED analysis^, the 
wrong bulk structure of Ref. |29| was used ia the scat- 
tering simulations. A GIXD measurements predicted 
also an inward relaxation of the topmost O-layer by 
—0.33 A and an outward relaxation of the underlying 
Zn-plane by +0.08 A. The occupancy probabilities were 
fitted, resulting in 1.3(!) and 0.7 for the first bilayer O 
and Zn, respectively. After considerably improved sam- 
ple preparation was achieved, the same aiithors reinves- 
tigated the O-terminated polar surface.Ej Best agree- 
ment with their GIXD data was now found for a struc- 
tural model where both, the upper O and Zn planes 
relax inwards by -0.19+0.02 A and -0.07+0.01 A, re- 
spectively, with occupancy factors of 1.0 in the oxygen 
plane and 0.75+0.03 in the underlying Zn plane. The 
inward relaxation of the O-layer has been confirmed by 
another surface X-ray diffraction measurementcj where 
Adi2=-0.24±0.06A and Ad 23 =+0.04±0.05 A was ob- 

Ab-initio calculations on polar slabsHEBB predict con- 
sistently for both surface terminations contractions for 
the first Zn-0 double-layer distance, with a larger in- 
ward relaxation at the O-terminated surface. 

In view of the above discussed discrepancies between 
different experimental and theoretical investigations, it is 
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our aim to provide a consistent set of fully converged cal- 
culations for the four main ZnO surfaces. We attempt to 
overcome the restrictions of previous theoretical studies 
such that the current study can be regarded as a refer- 
ence for perfectly ordered, defect-free surfaces. An ac- 
curate set of uniform theoretical data may then allow to 
discuss the differences between theory and experiment in 
terms of deviations between the model of ideal, unrecon- 
structed surfaces as assumed in the ab-initio simulations 
and the structure of the surfaces occurring in nature. In 
particular, for the polar surfaces this may give new in- 
sight into how these surfaces are stabilized. 



II. THEORETICAL DETAILS 
A. Method of calculation and bulk properties 

We have carried out self-consistent total-energy calcu- 
lations within the framework of density-functional the- 
ory (DFT)H3 The exchange and correlation effects were 
treatecLjplhin both, the local-density approximation 
(LDA)J£3S and the generalized-gradient approximation 
(GGA) where jsie used the functional of Perdew, Becke, 
and Ernzerhofi (PBE). 

Two different pseudopotential schemes were applied: 
For the study of the nonpolar surfaces we-used pseudopo- 
tentials of the Vanderbilt ultrasoft typeJUJ The electronic 
wave functions were expanded in a plane wave basis set 
including plane waves up to a cut-off energy of 25 Ry 
A conjugate gradient technique as described in Ref. |d] 
was employed to minimize the Kohn-Sham total energy 
functional. 

For the calculations on the polar surfaces we 
used normconserving pseudopotentialscia together with 
a mixed-basis consisting of plane waves and non- 
overlappipg localized orbitals for the 0-2p and the Zn-3d 
electronsJSI A plane- wave cut-off energy of 20 Ry was suf- 
ficient to get well converged results. To improve conver- 
gence in the presence of partly occupied bands, a Gaus- 
sian broadeningcfl with a smearing parameter of 0.1 eV 
was included. For several configurations representing 
nonpolar surfaces we repeated the calculations with the 
mixed-basis approach. No significant differences com- 
pared to the results from the ultrasoft-pseudopotential 
method could be seen. 

It is a well known shortcoming of LDA and GGA that 
both predict the Zn-<i bands to be roughly. 3 eV too 
high in energy as compared to experiment .ElfCj In conse- 
quence, the Zn-ci states hybridize stronger with the O—p 
valence bands, thereby shifting them unphysically close 
to the conduction band. The underestimate for the band 
gap is therefore even more severe in ZnO than in other 
semiconductors. In our calculations we obtained band 
gaps of 0.78 eV and 0.74 eV with LDA and PBE, respec- 
tively, as opposed to the experimental value of 3.4 eV. 
The band gap and the position of the Zn-<i bands can 
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FIG. 3. Schematic top and side view of the surface ge- 
ometry of the nonpolar (1120) surface. The same representa- 
tion as in Fig. ^ is used. The glide planes are indicated by 
dashed-dotted lines. 



be improved significantly, if a self-interaction correction 
(SIC) is used.E3 Usually SIC calculations are very de- 
manding, but if the, SIC effects are incorporated into 
the pseudopotentialuS, the additional calculational cost is 
modest. Unfortunately, the SIC pseudopotential scheme 
does not improve the structural properties of ZnOcS and 
also causes .some problems when accurate atomic forces 
are needed.EJ Therefore, and since we are mostly inter- 
ested in accurate relaxed geometries of the surfaces and 
not so much in their electronic structure, we have omitted 
the use of SIC in our calculations. 

The computed structural parameters for bulk ZnO 
are shown in Table H. Mixed-basis and ultrasoft- 
pseudopotential calculations give the same results within 
the accuracy displayed in Table Q. As is typical for the 
functionals, LDA underestimates the lattice constants by 
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(a) Top view: 




(b) Side view: 
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FIG. 4. Schematic top and side view of the polar 
Zn-terminated (0001) surface. The same representation as 
in Fig. ^ is used. 



1-2 %, and GGA overestimates them by roughly the same 
amount. The c/a-ratio strongly influences the internal 
parameter u. If u — 1/4 + a 2 /3c 2 , all nearest-neighbor 
bonds are equal. Since the c/a-ratio is slightly over- 
estimated in our calculations, we get u— values that are 
slightly smaller than observed in experiment. 

The construction of appropriate supercells for the 
study of the surfaces will be detailed in the following 
subsection. All atomic configurations were fully relaxed 
by minimizing the atomic forces using a variable-metric 
scheme. l£I Convergence was assumed when the forces on 
the ions were less than 0.005 eV/A. 



B. Surfaces, slab structures, and the stability 
problem 



faces a dipole correctionBB was used to prevent artificial 
electrostatic interactions between the repeated units. To 
simulate the underlying bulk structure, the slab lattice 
constant in the direction parallel to the surface was al- 
ways set equal to the theoretical equilibrium bulk value 
(see Table |). 

The nonpolar surfaces are obtained by cutting the crys- 
tal perpendicular to the hexagonal Zn- and O-layers (see 
Fig. 0). In both cases, for the (1010) and the (1120) 
planes, two equivalent surfaces are created so that always 
stoichiometric slabs with the same surface termination on 
top and on bottom can be formed. 

The (1010) surface geometry is sketched in Fig. |^. 
Each surface layer contains one ZnO dimer. The dimers 
form characteristic rows along the [1210] direction which 
are separated by trenches. Slabs with 4-20 atomic lay- 
ers were used, thus containing up to 40 atoms, and 
the Brillouin-zone of the .Sitpercell was sampled with a 
(4x2x2) Monkhorst-PackEil k-point grid. No differences 
were found when going to a (6x4x2) mesh. 

The surface layers of the (1120) surface are built up 
by two ZnO dimers which form zig-zag lines along the 
surface (see Fig. ^). The two dimers are equivalent and 
are related by a glide plane symmetry. This symmetry k, 
not destroyed by the atomic relaxations of the surfaced 
The slabs in our calculations were built of 4-8 atomic 
layers_with up to 32 atoms, and a (2x2x2) Monkhorst- 
Packcil k-point mesh was used. Again, a denser (4x4x2) 
mesh did not alter the results. 

Cleaving the crystal perpendicular to the c-axis (see 
Fig. |l|) always creates simultaneously a Zn- and an O- 
terminated polar (0001) and (0001) surface, respectively. 
If we only consider cuts where the surface atoms stay 
3-fold coordinated, all slabs representing polar surfaces 
are automatically stoichiometric and are inevitably Zn- 
tcrminated on one side and O-terminated on the other 
side. Figure |J sketches the characteristic sequence of Zn- 
O double-layers of the polar slabs. In our calculations 
slabs with 4-20 Zn-0 double-layers were used, thus con- 
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All surfaces were represented by periodically repeated 
slabs consisting of several atomic layers and separated 
by a vacuum region of 9.4 to 12.4 A. For the polar sur- 



FIG. 5. Schematic illustration of the stacking sequence of 
the polar slabs. A charge transfer of 5 = (1— 2m) c ps 2/4 has 
to occur to stabilize the polar surfaces. 
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FIG. 6. Schematic illustration of the band structure after 
electrons have moved from the O- to the Zn-terminated sur- 
face of the slab. Depending on the band gap and the thickness 
D of the slab, a residual electric field remains inside the slab. 

taining 8-40 atoms, k-point corpergence was achieved 
with a (6x6x1) Monkhorst-PackEil grid, and tests with 
up to (12x12x1) k-points were made. 

Each Zn-0 double-layer in Fig. |I] exhibits a dipole 
moment perpendicular to the surface. If we assume for 
simplicity a purely ionic model for ZnO and assign the 
fixed formal charges -\-Ze and — Ze to the Zn- and O- 
ions, respectively, then a slab of N double-layers will 
exhibit a dipole moment of m = N Ze (1 — 2u) c/2 (see 
Fig. ||). This corresponds to a spontaneous polarization 
of P s = Ze (1 — 2u) which is independent of the thickness 
of the slab. If the external electric field is zero, inside 
the slab an electric field of E = —4irP s will be present. 
Therefore, no matter how thick we choose our slab, the 
inner part will never become bulk-like, and the surface 
energy, defined as the difference between the energy of 
the slab and the energy of the same number of atoms in. 
the bulk environment, will diverge with slab thickness SB 
Thus, the polar surfaces are not stable. 

On the other hand, it can easily be seen that if we 
modify the charge in the top and bottom layer of the slab 
from ±Ze to ±(Z-S)e with 5=(l-2u)Z « Z/4, then 
the dipole moment of the slab will become independent of 
the slab thickness and the internal electric field vanishes. 
This charge transfer is equivalent to applying an external 
dipole which compensates the internal electric field. 

For most polar surfaces the rearrangement of the 
charges is accomplished by a modification of the surface 
layer composition with respect to the bulk. If this does 
not occur, the internal electric field will 'tilt' the band 
structure by which the upper edge of the valence band 
close to the O-terminated surface will become higher in 
energy than the lower edge of the conduction band at 
the Zn-terminated face (see Fig. The slab can now 
lower its energy (thereby reducing the internal electric 
field) by transferring electrons from the valence band at 
the O-terminated side to the conduction band at the Zn- 
terminated face. This will happen 'automatically' in any 
self-consistent electronic structure calculation that makes 
use of a slab geometry. This is what is usually referred 
to as 'the metallization of polar surfaces'. 



However, one problem still remains: electrons move 
from the O- to the Zn-terminated surface until the upper 
valence band edge at the O-terminated side has reached 
the same energy as the lower edge of the conduction band 
at the Zn-terminated face as sketched in Fig. |^. In this 
situation, the internal electric field is not fully removed 
for a finite slab with thickness D. The residual electric 
field depends on the band gap and vanishes only with 
l/D. In our calculations we found that for slabs with up 
to 6 Zn-0 double-layers the residual electric field is still 
so strong that the slabs are not stable. There is no energy 
barrier when the O and Zn layers arc shifted simulta- 
neously and rigidly toward each other. Therefore, to get 
good converged results for the surface geometries and en- 
ergies very thick slabs have to be used which makes the 
investigation of the polar surfaces computationally very 
demanding. Ideally, one should calculate all quantities 
of interest for different slab thicknesses D and extrapo- 
late the results to l/D — ► 0. In the present study we 
obtained the relaxations of the surface layers (see Fig. ||) 
as well as the cleavage energy of the polar surfaces (see 
Fig. ^) by extrapolating the results of slab calculations 
containing up to 20 Zn-0 double-layers. 

III. RESULTS AND DISCUSSION 
A. The nonpolar (1010) and (1120) surfaces 

The nonpolar wurtzite (1010) surface and the closely 
related zincblendc (110) surface have been studied exper- 
imentally and theoretically for a wide range of III-V and 
II- VI compound semiconductors. It was found that all 
surfaces show the same basic relaxations with the sur- 
face cation moving inwards and the anion staying above, 
resulting in a tilt of the surface anion-cation dimers, 
and the magnitude of the relaxation is determined by a 
competitionJaetween dehybridizaton and charge transfer 
eSectsJBWB 

At the surface (this applies also to the (1120) surface), 
the coordination of the surface atoms is reduced from 
4-fold to 3-fold, thereby creating an occupied dangling 
bond at the anion and an empty dangling bond at the 
cation. Two limiting cases may now be distinguished: In 
a dominantly covalent bonded compound the cation will 
rehybridize from sp 3 to sp 2 and will move downwards 
until it lays nearly in the plane of its three anion neigh- 
bors. The anion stays behind (often even an outward 
relaxation is observed) tending toward p-like bonds to 
its neighbors. The result is a strong tilt of the surface 
anion-cation dimer (up to 30° are observed) with only a 
small change of the bondlength. In a dominantly ionic 
solid, electrostatics prevails over dehybridization effects. 
To obtain a better screening, both, anion and cation, 
move toward the bulk. The tilt of the anion-cation dimcr 
will be small but the bond length can be significantly re- 
duced. Therefore, the relaxation of the surface dimers 
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(1010) surface 
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10X 





LDA 


PBE 


LDA 




PBE 


Al,± 


+0.106 a 


+0.100 a 


+0.076 a 




+0.073 a 


A 2 ,± 


-0.041 a 


-0.038 a 


-0.016 a 




-0.015 a 


Bulk 


0.0 






0.0 




Ai, H 


0.6531 c 


0.6539 c 


0.6506 c 




0.6516 c 


A 2 , y 


0.6243 c 


0.6231 c 


0.6230 c 




0.6221 c 


Bulk 


(1-u 


)c 


(1 


-u) 


c 


Ai, x 






0.083 a 




0.077 a 








0.020 a 




0.016 a 


Bulk 








0.0 




di2,± 


0.1445 a 


0.1447 a 


0.4093 a 




0.4089 a 


Bulk 


0.6328 a 


0.6337 a 
fa 


0.5215 a 


a/2 


0.5222 a 


d\2,y 


0.5355 c 


0.5357 c 


0.5259 c 




0.5266 c 


d23,y 


0.5013 c 


0.5017c 


0.5014 c 




0.5009 c 


Bulk 


c/2 






c/2 










0.4381 a 




0.4399 a 


d2,x 

Bulk 






0.5515 a 


fa 


0.5556 a 



TABLE II. Summary of the structural relaxations of the 
first two surface layers for the nonpolar (1010) and (1120) 
surfaces. The definitions of the independent structural pa- 
rameters are shown in Fig. and ^. All relaxations are given 
in fractions of the theoretical bulk lattice constants a and c 
(see Table |j|) . The rows labeled 'Bulk' are the corresponding 
values for the unrelaxed surface. 







C b,\\ 


C* B (Zn) 


C B (0) 


LDA, this study 


10.7° 


-6.7 


-2.8 


-3.2 


PBE, this study 


10.1° 


-7.2 


-3.1 


-3.4 


LEED, Ref. | 


12°±5° 


-3+6 






LDA+pw, Ref. |(] 


11.7° 


-6.0 






LDA+Gauss, Ref. [i]] 


3.6° 


-7.9 


-5.2 


-2.7 


HF, Ref. [l| 


2.3° 


-7.2 


-3.6 


-3.4 


B3LYP, Ref. [| + || 


5.2° 


-4.9 


-2.9 


-0.5 



TABLE III. Tilt angle lu of the surface dimer (see Fig. |) 
and relative bond length contraction Cb (in % of the cor- 
responding bulk value) of the surface bonds for the nonpolar 
(1010) surface in comparison with LEED experiment and pre- 
vious calculations. Cb,|| refers to the Zn-O dimer bond paral- 
lel to the surface, Ca(Zn) to the back bond of zinc to oxygen 
in the second layer, and Cb(O) to the respective back bond of 
the surface O atom. Bulk values of the surface and the back 
bonds are uc and ((1/2 — u) 2 + a 2 /3c 2 ) c, respectively. 




1 23456789 10 
surface layer 



bond length 
contraction C„ 



123456789 10 
surface layer 



FIG. 7. Deep-layer relaxations for the nonpolar (1010) 
surface calculated with a 20 layer slab. Plotted are the tilt 
angle lu (in degrees) and the in-plane bond length contraction 
C B | (in %) of the Zn-0 dimers as a function of the distance 
from the surface. Only the PBE results are shown, the results 
from the LDA calculations are essentially identical. 



directly reflects the covalency or ionicity of the chemical 
bond in the compound of consideration. 

Our results for the relaxation of the (1010) surface are 
given in Tables |l| and [II. All lengths are expressed as 
fractions of the theoretical lattice parameters given in 
Table ffl. Using these dimensionless relative quantities no 
significant differences between the LDA and GGA calcu- 
lations can be seen. For two structural parameters the 
decay of the surface relaxations into the bulk is illustrated 
in Fig. 0. Compared to the topmost surface layer, the tilt 
angle lu and the in- plane bond length contraction Cbj| of 
the Zn-0 dimers are already much smaller in the second 
and the subsequent layers, but still significant deviations 
from the bulk structure can be seen as deep as five or six 
layers below the surface. 

The relatively small angle of lu s=s 10° for the tilt of the 
surface Zn-0 dimer together with the Zn-0 bond con- 
traction of Cbji ~ 7 % confirms that the chemical bond in 
ZnO is highly ionic but with significant covalent contribu- 
tions. A tilt of 10° is at the lower boundary of what has 
been observed for other III-V and II- VI compounds.^ 
Only the nitride-semiconductors show tilt angles that are 
similarly small.E3 

The calculated surface relaxations in Table || and [II 
agree very well with the DFT-LDA study of Ref. ^ and 
with the results from the LEED analysis B Relative to 
the central layer of the slab we find a downward re- 
laxation of the surface atoms of Ari^(Zn)=— 0.36 A and 
Ac?_l(0)=— 0.04 A with a shift parallel to the surface of 
Ad||(Zn)=0.18A compared to Ad ± (Zn)=-0.45±0.1 A, 
Ari_L(O)=-0.05±0.lA, and Ad||(Zn)=0.1±0.2A from 
the LEED experiment. 

Rotation angles of lu—2°-5° seem anomalously small 
in the context of what is known for other compounds. 
Even for the vpy ionic A1N a tilt angle of lu=7.5° has 
been reported.E2l The smaller relaxations obtained in 
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0.12 





w 


Cb,|| 


C B (Zn) 


Cb(O) 


LDA, this study 


7.6° 


-5.8 


-1.4 


-1.7 


PBE, this study 


7.4° 


-6.4 


-1.5 


-1.8 



TABLE IV. Tilt angle ui of the surface dimer (see Fig. |) 
and relative bond length contraction Cb of the surface bonds 
for the nonpolar (1120) surface. The same notation as in 
Table III is used. 



Rcf. |17|-[19| may be due to not fully converged calcula- 
tions. Very thin slabs were partly used or only the first 
one or two surface layers were relaxed. In Ref. [l^ no 
relaxation of the Zn ions parallel to the surface was al- 
lowed. Also the convergence of the localized basis sets 
employed in these studies and the k-point sampling may 
have been a problem. As a test we did a slab calculation 
where we fixed the positions of the atoms at the bulk po- 
sitions and allowed only the first surface layer to relax. 
The tilt angle u) then reduces to roughly half of its value. 
Also coarsening the k-point mesh results in changes of 
2°-3° in u>. Since we did our calculations with two very 
different pseudopotential approaches we can exclude any 
bias caused by the use of pseudopotentials. 

Table H and [fv| show our results for the relaxation of 
the (1120) surface. The atomic displacements are of the 
same order of magnitude as has been found for the (1010) 
surface. Again, no significant differences between LDA 
and GGA calculation can be seen. The tilt of the surface 
dimers of 7.5° and the reduction of the Zn-0 dimer bond 
length of about 6 % fits nicely into the picture of ZnO 
being at the borderline between ionic and covalent solids. 

In a hybrid B3LYP studjO much smaller relaxations 
for the (1120) surface were reported. However, in this 
study only three degrees of freedom per surface layer were 
relaxed. The authors claimed that the position of the 
Zn and O ions are constrained by symmetry. This is 
not correct. From the two Zn-0 dimers in each surface 
layer, the atoms of one dimer can move freely in all three 
Cartesian directions, leading to six degrees of freedom 
per surface layer (see Fig. ||) ■ The position of the second 
dimer is then determined by the glide plane symmetry 
(see also Ref. p3|). 



B. The polar (OOOl)-Zn and (OOOl)-O surfaces 

In Figure || we have plotted the calculated distances 
between the topmost surface layers of the polar (0001) 
and (0001) surfaces as a function of the slab thickness 
D. As expected from the thickness dependence of the 
residual electric field inside the slab, the 1/D plots re- 
veal a nice linear behavior for the interlayer distances. 
By extrapolating 1/D — > 0, all distances may now be 
obtained in the limit of a vanishing internal electric field. 

The extrapolated results for the relaxations of the po- 



0.08 



0.06 



0.42 




2016 12 10 8 
double-layers N 



0.38 



2016 12 10 8 
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FIG. 8. First three interlayer distances (see Fig. ^) for 
the polar Zn-terminated (0001) and the oxygen terminated 
(0001) surface calculated with different slabs containing N 
Zn-O double-layers and using the PBE functional. All dis- 
tances are given in fractions of the theoretical bulk lattice 
constant c (see Table |l]) and are plotted vs. 1/N. The extrap- 
olation 1 /N — > gives the surface relaxations for a vanishing 
internal electric field. 

lar surfaces are summarized in Table |v| and [v|. Very 
good agreement with the results of previous ab-initio cal- 
culations is found. In general, all double-layers are con- 
tracted and the distances between the double-layers are 
increased relative to the bulk spacings. For finite slabs, 
the residual internal electric field further amplifies this 
characteristic relaxation pattern. 

The largest relaxation is found for the O-terminated 
surface where the outermost double-layer distance is com- 
pressed by «50%. This agrees reasonably well with the 
results of the X-ray experimentsc3'lEiL3 where a contrac- 
tion of 40%, 54%, and 20%. were found. On the other 
hand, from LEED analysisES and LEISl^I measurements 
it was concluded that the Zn-0 double-layer spacing for 
the O-terminated surface ispclose to its bulk value. The 
recent finding of W611 et al.cl may perhaps help to solve 
this contradiction. With helium scattering it was shown 
that after commonly used preparation procedures the O- 
terminated surfaces are usually hydrogen covered. To 
test how much hydrogen may influence the surface re- 
laxations, we repeated a calculation where we adsorbed 
hydrogen on top of the O-terminated side of the slab. We 
find that in this case the outermost Zn-0 double-layer 
expands again, and the Zn-0 separation goes back close 
to the bulk distance. A similar result was also reported 
by Wander and Harrison 

For the Zn-terminated surface there is a clear dis- 
crepancy between theory and experiment. All calcula- 
tions predict consistently a contraction of the first Zn- 
O double-lajjer of 20-30%, whereas in experiment no 
contractionEd or even an mi.tward relaxation of the top- 
most Zn-layer is found.Ejo This may indicate that the 
'metallization' used in all theoretical studies is not the 
adequate model to describe the polar JZji-terminated sur- 
face. Recently Dulub and DieboldEll proposed a new 
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4.0 



(OOOl)-Zn surface 
LDA PBE 



(OOOl)-O surface 
LDA PBE 



di2 
d 2 3 

d,34 
'hr, 
ds6 

Bulk 



0.0952 
0.3947 
0.1172 
0.3811 
0.1187 



0.0883 
0.3995 
0.1132 
0.3857 
0.1186 



0.0594 
0.4022 
0.1044 
0.3817 
0.1194 



0.0645 
0.3962 
0.1107 
0.3784 
0.1251 



(A — u) C I UC 



(h—u)c I UC 



TABLE V. Summary of the surface relaxations for the po- 
lar Zn-terminated (0001) and the O-terminated (0001) sur- 
face (see Fig. ^|). All distances are given in fractions of the 
theoretical bulk lattice constant c (see Table ||). 







(0001)-Zn surface 


(0001)-O surface 






Adi2 


Ad 23 


Adi2 


Ad 23 


LDA, this s 


tudy 


-22% 


+5.1% 


-51% 


+4.7% 


PBE, this s 


tudy 


-27% 


+5.3% 


-47% 


+4.5 % 


B3LYP, Ref. §| 


-23% 


+3.5% 


-41% 


+3.0 % 


GGA, Ref. 


27] 


-31% 


+7.0% 


-52% 


+6.5 % 


GGA, Ref. 




-25% 




-41% 





TABLE VI. Relaxation of the surface layers of the polar 
ZnO surfaces in comparison with previous ab-initio calcula- 
tions. 



stabilization mechanism for the Zn-terminated surface. 
With scanning tunneling microscopy (STM) they found 
that many small islands with a height of one double-layer 
and many pits one double-layer deep are present on the 
(0001)-Zn surface. Assuming that the steps edges are 
O-terminated, an analysis of the island and pit size dis- 
tribution yielded a decrease of surface Zn concentration 
of roughly 25%. Such a reduction of Zn atoms at the 
surface would be enough to accomplish the charge trans- 
fer needed to stabilize the polar surface. It would not 
be in contradiction with the observed p(l x 1) LEED pat- 
tern since a long range correlation between the different 
terraces remains. The missing of 25 % of the Zn atoms 
was also obtained by JedrecjH as best fit of their GIXD 
data. 

A structure where the surface is stabilized by many 
small islands and pits with a Zn deficiency at the step 
edges is, of course, far away from the model of a clean, 
perfectly ordered (0001)-Zn surface used in the theoret- 
ical calculations. Basically all surface Zn-atoms will be 
next to a step edge, and therefore very different relax- 
ations may occur. Unfortunately, it is presently out of 
the reach of our ab-initio method to do calculations on 
slabs representing such an island and pit structure. 

For the O-terminated surface, on the other hand, 



2.8 
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FIG. 9. Similar plot as Fig. ^| for the cleavage energy of 
the polar ZnO surfaces. Shown are the results of the PBE 
calculations. 



the STM measurements show a very different picture. 
Smooth and flat terraces separated mostly by a two 
double-layer step are observed. The number of single 
double-layer steps was by far not large enough to ac- 
count for a similar stabilization mechanism as for the 
Zn-terminated surface. 



C. Surface/cleavage energies 

For the nonpolar surfaces we can obtain directly the 
surface energy from our slab calculations since the slabs 
are always terminated by the same surface on both sides. 
This is not possible for the polar surface since inevitably 
both surface terminations are present in a slab calcula- 
tion. Only the cleavage energy of the crystal is well de- 
fined. To be able to compare the relative stability of the 
nonpolar and polar surfaces, we will discuss in the fol- 
lowing only the cleavage energies. The surface energies 
of the nonpolar surfaces are just given by half of their 
cleavage energy. 

Like the interlayer distances, the 1/D-plot of the cleav- 
age energy for the polar surfaces in Fig. ^ exhibits a sim- 
ple linear behavior. As can be seen, the cleavage energy 
does not change too much with the slab thickness so that 
moderate slab sizes would be sufficient to obtain reason- 
able converged results. 

The extrapolated values for the cleavage energy of the 
polar surfaces together with the results for the nonpolar 
faces and the findings of previous studies are summarized 
in Table VII. The nonpolar (1010) surface is the most 
stable face of ZnO with the lowest cleavage energy. But 
the energy of the (1120) surface is only slightly higher. 
The cleavage energy for the polar surface is roughly a 
factor of two larger than for the nonpolar surfaces. This 
is surprisingly low compared to what has been found 
in other systems, for example MgO, where a 'metalliza- 
tion' was also assumed as stabilization mechanism for 
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Ert 



(1010) surface: 

LDA, this study 2.3 0.37 

PBE, this study 1.6 0.37 

LDA+pw, Refjio] 1.7 0.37 

B3LYP, Ref. [jj + || 2.3 

HF, Ref. |l| 2.7 0.38 

Shell model, Ref. |l| 2.0 

(1120) surface: 

LDA, this study 2.5 0.29 

PBE, this study 1.7 0.30 

(0001)/(000l) surface: 

LDA, this study 4.3 0.28 

PBE, this study 3.4 0.28 

B3LYP, Ref. |^ 4.0 

TABLE VII. Cleavage energy i? c i eav (in J/m 2 ) and relax- 
ation energy i? re iax (in eV per surface Zn-O dimer) for the 
different ZnO surfaces and in comparison with previous cal- 
culations. 



the polar surfacesEa Therefore, for ZnO the 'metalliza- 
tion' mechanism can well compete with other stabiliza- 
tion mechanisms like reconstructions or randomly dis- 
tributed vacancies and can not be ruled out by energetic 
considerations alone. 



The LDA and GGA results in Table VII show the same 



ordering for the cleavage energies of the different surfaces 
but the absolute GGA energies are roughly 30 % lower 
than the LDA results. This is a well known improve- 
ment of the GGA, where a much better description of 
the rapidly decaying charge density into the vacuum re- 
gion is achieved. The cleavage energies agree well with 
previous theoretical results as given in the Table. Sur- 
prisingly, the results of the hybrid B3LYP studies are 
much closer to LDA than to the GGA results. Interest- 
ingly, the relaxation energy is roughly the same for all 
surfaces when normalized to one Zn-0 pair. This means 
that despite the partially filled bands at the polar sur- 
faces, the strength of the relaxation is almost the same 
as for the isolating nonpolar faces. 



IV. SUMMARY AND CONCLUSIONS 

A first-principles density-functional pseudopotential 
approach was used to determine the fully relaxed atomic 
structures and the surface/cleavage energies of the non- 
polar (1010) and (1120) surfaces and the polar Zn- 
terminated (0001) and the O-terminated (0001) basal 
surface planes of ZnO. 



The main results of the presented investigation are an 
extensive set of reliable data for the structural parame- 
ters and the energetics of the various ZnO surfaces within 
the LDA and the PBE approximation, which we consider 
to be a reference for future studie s (se e in particular the 
compilations in Tables |l[ |v|, and VII ) . 

For the nonpolar surfaces we could resolve the discrep- 
ancy between experiment and several previous ab-initio 
studies by showing that if calculations are carefully con- 
verged a moderate tilt of the Zn-0 surface dimers with a 
relatively strong contraction of the dimer bond length is 
obtained. Such a relaxation pattern is typical for rather 
ionic compounds but with strong covalent contribution 
to the chemical bonding. Our results are in line with 
LEED analysis and fit very well the systematic trends 
that are observed for other more or less ionic II- VI and 
III-V semiconductors. 

The polar surfaces can only be stable if a rearrange- 
ment of charges between the Zn- and the O-terminated 
surfaces takes place. In our calculations the polar sur- 
faces were stabilized by allowing the electrons to move 
from the (0001)-O to the (0001)-Zn surface, thereby 
quenching the internal electric field. Nevertheless, even 
for thick slabs a finite residual electric field is present in- 
side the slabs, which affects the results for the structural 
parameters and the surface energy. To get well converged 
results in the limit of a vanishing internal electric, we re- 
peated all calculations with slabs consisting of different 
numbers of Zn-0 double-layers and extrapolated the re- 
sults to the limit of an infinite thick slab. 

For both polar surfaces we obtain a strong contraction 
of the outermost double-layer spacing. This agrees well 
with experiments for the O-terminated surface but not 
for the Zn termination, indicating that the electron trans- 
fer may be not the adequate model to describe the stabi- 
lization mechanism of the polar Zn-terminated surface. 
Since this is consistently predicted by all calculations, it 
is very likely that other mechanisms, such as defect for- 
mation, hydroxylation and/or the mechanism proposed 
by Dulub and Diebold might stabilize the (0001)-Zn sur- 
face. 

Concerning the surface energies, we find very similar 
values for the two nonpolar surfaces with a slightly lower 
value for the (1010) surface. The cleavage energy for the 
polar surfaces is predicted to be roughly a factor of two 
larger than for the (1010) face. 
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